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We present a family of discrete breathers, which exists in a nonhnear polarizabihty model of 
ferroelectric materials. The core-shell model is set up in its non-dimensionahzed Hamiltonian form 
and its linear spectrum is examined. Subsequently, seeking localized solutions in the gap of the 
linear spectrum, we estabUsh that numerically exact and potentially stable discrete breathers exist 
for a wide range of frequencies therein. 

In addition, we present nonhnear normal mode, extended spatial profile solutions from which 
the breathers bifurcate, as well as other associated phenomena such as the formation of phantom 
breathers within the model. The full bifurcation picture of the emergence and disappearance of the 
breathers is complemented by direct numerical simulations of their dynamical instability, when the 
latter arises. 



INTRODUCTION 

The study of discrete breathers has gained con- 
siderable momentum over the past decade. This has 
been motivated chiefly by their experimental obser- 
vation in a wide variety of physical systems. Such 
exponentially localized in space and periodic in time 
nonlinear waveforms have been identified in complex 
electronic materials such as halide-bridged transition 
metal complexes as, e.g., in they have also been 
proposed to be relevant in the formation of denatura- 
tion bubbles in the DNA double strand dynamics as 
discussed, e.g., in Such localization has also been 
discussed in the context of ultracold Bose-Einstein 
condensates immersed in optical lattices [3j. Addi- 
tional prototypical frameworks for the experimental 
formation of such modes include the nonlinear op- 
tics of evanescently coupled waveguide arrays, and 
biased photorefractive crystals; a recent review can 
be found in [2] . Localized modes have also been ob- 
served in micromechanical cantilever arrays [5j], in- 
coupled torsion pendula 0] , in electrical lattices and 
transmission lines f^, in granular crystals of beads 
interacting through Hertzian contacts ^ , and in lay- 
ered antiferromagnetic saniples such as those of a 
(C2H5NH3)2CuCl4 crystal 

In the context of ferroelectrics, a prototypical 
model for their nonlinear dynamical description can 
be traced back to the work of Bilz et al. Ill |. 



roelectrics including the temperature dependence of 
soft-mode frequencies and dielectric constants [121] . 
The theoretical analysis of this model has generally 
required considerable simplifications and assump- 
tions [llj. Yet, given the fundamental ingredients of 
nonlinearity and discreteness and their interplay rel- 
evant for such systems, discrete breathers (or intrin- 
sic localized modes) have been discussed as poten- 
tially arising within this model. They are associated 
with the formation of polar nanoregions indicating a 
crossover between and coexistence of displacive and 
relaxor behavior of ferroelectrics 
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It was subsequently shown that, independently of 
structural differences, this nonlinear polarizability 
model is applicable in general to displacive-type fer- 



In this paper, our aim is to undertake a system- 
atic numerical study of the existence and stability 
properties, as well as the nonlinear dynamics, of 
breathers within this polarization model for a pro- 
totypical set of parameter values, relevant to re- 
alistic contexts (ferroelectric materials). We illus- 
trate that breathers not only exist but can be lin- 
early stable within the dispersion gap of the under- 
lying linear system. In the process, we also explore 
some of the interesting symmetry-breaking (pitch- 
fork) bifurcations that give rise to the emergence of 
these waveforms within the nonlinear polarizability 
model and some of the (second harmonic genera- 
tion) resonance phenomena that they incur. Our 
presentation is structured as follows: In section II, 
we give the theoretical description of the model and 
its underlying linear spectrum. In section III, we 
present our numerical results for the existence, sta- 
bility and dynamical properties of the obtained dis- 
crete breathers. Finally, in section IV, we summa- 




FIG. 1. The schematic of the pseudo-one dimensional 
nonlinear polarizability model for displacive-type ferro- 



electrics. A diatomic lattice with core {ui^ 



U2- 



and 



shell (wij, V2,j) displacements along with the connecting 
springs is shown for two neighboring sites i and i + 1. 



rize our findings and present some conclusions and 
suggestions for future studies. 



THEORETICAL MODEL 

The nonlinear polarizability model describes a fer- 
roelectric material by grouping atoms in the full sys- 
tem into a "pseudo-NaCl" structure, and then pro- 
jecting it onto the direction of ferroelectric polariza- 
tion ll| , simplifying a 3-D lattice to a diatomic I-D 
lattice that retains many of the features of the orig- 
inal model. Every lattice site has two atoms, each 
with a core and a shell (Fig. 1). The displacements 
from equilibrium of cores one and two at lattice site 
i are ui.i and U2,i, and the displacements of the cor- 
responding shells are vi^i and V2,i- The Hamiltonian 
of the model is given by 
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TABLE I. Physical parameter values in Eqs. (2) and 
(3). 



Here mi, Wei, W2 and 111^2 denote the core and 
shell masses of atoms 1 and 2, respectively. The 
shell-shell, core-core, and internal core-shell har- 
monic force constants are given by /, and 52*^ 
for i = 1,2. The anharmonic force constant 54 takes 
into account the nonlinear polarizability of atom 1. 
Typical values of the physical parameters,_which 
were taken to match those of KTaOs in [12[ are 
shown in Table HI 

It is convenient to rescale the Hamiltonian ([IJ to 
the dimensionless form, so that the motion of shell- 1 
is of order unity. This can be achieved by making 
the substitutions t — tcT, u = XclJ^, and v = Xc'y, 
and H = ecH, where 
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H = T + V, 



where 
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Now, for convenience, we drop the hats and rela- 
bel /i and 7 as u and w; we will make sure, from 
this point forward, to specify when we are referring 
to the dimensional form. Note here that any so- 
lution obtained in the dimensionless system can be 
re-dimensionalized by multiplying length by Xc, time 
by tc, or energy by Cc] these are shown in Table HI 
This gives us the resulting non-dimensional Hamil- 
tonian 



H = T + V, 



(5) 
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where 

T = i ^{rhiul i + m2u\, + vf^, + me2V2,i)^ (6) 

i 
i 

+ D{v2.i - Vi,if + D{v2,i - Vl.i+lT 

+ \{vi,^ - u^,i)% (7) 
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mei mei 
This results in the fohowing equations of motion 

miMi,i = Ci(ui,i_i + ui,i+i - 2ui,i) 

- {vi,i - Ui,i) + {vi,i - Ui,i)^, 
m2U2,i C2{U2,1-1 + U2,i+l - 2u2,i) 

+ a2{v2.i - U2,i), 
Vi^i = D{v2,n-1 + V2,n " 2wi^„) 

+ ivi,n - Ui,n) - (Vi^n - ■til,„)^, 
me2V2,i = D{vi^n + f l,n+l " 2^2, n) 

- a2{v2,n ~ U2,n)- 

The schematic of this model is shown in Fig. [TJ 
Given the effective four-component nature of the 
model and the computational expense of perform- 
ing the (up to a prescribed accuracy) exact nu- 
merical computation amounting to the identifica- 
tion of the breather solutions discussed below, we 
have restricted our consideration to finite size lat- 
tices with N sites and periodic boundary conditions 
{ujfi = Uj^N+i, Vjfi = Vj^N+i, j = 1,2). Typically 

= 20 has been used in the results shown below, 
although we have ensured that this choice does not 
significantly affect the nature of our existence and 
stability results. 

Note also that this system is invariant under trans- 
lations, so we have chosen (for ease of visualization) 
to shift each solution so that the spatial average 
of core-1 is zero. Lastly, this system is clearly in- 
variant under parity changes, so if the vector x = 
(ui^n, U2,„, wi,n, W2,„) reprcscuts a solution, then so 
does —X. 
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FIG. 2. The linear spectrum of both linearizations (i.e. 
stable and unstable equilibria) with four bands and three 
frequency gaps. Note that the bottom two bands over- 
lap. The right two panels show the eigenmodes from the 
bottom edge of the top band, and the top edge of the 
bottom band. In these panels, core/shell-l are shown 
in blue circles, while core/shell-2 are shown as red di- 
amonds. Cores are the filled markers, while shells are 
open markers. 



Linear Spectrum 

The above system possesses both stable and un- 
stable uniform state equilibria. Of the former kind 
are the solutions with ui^i = 0, U2.i = Vj^i = ±1, 
while of the latter kind are the vanishing solutions 
with Uj^i = Vj_i = 0. Equations ([9]) were linearized 
about the equilibrium solutions to obtain the lin- 
earized equations x — J{x*)x, where J{x) is the 
Jacobian matrix of the system ([3]), and x* is the 
equilibrium vector being linearized around. The lin- 
ear spectrum changes depending on the particular 
form of X about which we linearize. 

More specifically, we looked for solutions of the 
form uin = exp(iwi), U2,„ = ^2,™ exp(ia;t), 
vi.7i — exp(za;t), and U2,n — J^2,n exp(iwt), and 
then solved the ensuing eigenvalue problem. The lin- 
ear spectrum of both the stable and unstable equi- 
libria is shown in Fig. [21 There are four bands and 
three frequency gaps (the bottom two bands over- 
lap). The first band has frequencies in the range 
of - 5.894 X 10"'^, and the dominant sublattice 
is that of corc-1. The second band has frequencies 
in the range of 5.438 x 10"^ - 11.57 x 10"^ and 
the dominant sublattice is that of core-2. The third 
band has frequencies in the range of 0.2589 - 0.5473 
Hz, and the dominant sublattice is that of shell- 1. 
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Lastly, the fourth band has frequencies in the range 
of 1.1949 - 1.2885 Hz, and the dominant sublattice 
is that of shell- 2. 



NUMERICAL EXISTENCE, STABILITY AND 
NONLINEAR DYNAMICS OF DISCRETE 
BREATHERS 

We use a Newton-Raphson solver to identify the 
zeros of a Poincare map for periodic orbits of the dis- 
crete breather type; for a relevant review of compu- 
tational methods for the identification of such states, 
see, e.g., [3]. Namely, we solve for the zeros of the 
map P{xo) = xq — x{T, xq), where x(T, xq) is a solu- 
tion satisfying the equations of motion with ini- 
tial condition xq, evaluated at time T. Subsequently, 
a numerical continuation was performed using the 
breather frequency as the continuation parameter. 
At each continuation step we calculated the stabil- 
ity of the periodic solution in the form of its Floquet 
multipliers by computing the eigenvalues of the mon- 
odromy matrix V{T), which is obtained by integrat- 
ing the variational equation V — J{x{t))V along 
the periodic solution x for one full period, with an 
initial condition of the identity matrix. Floquet mul- 
tipliers with a magnitude greater than one are asso- 
ciated with an unstable solution in our Hamiltonian 
system (where if A is a multiplier, then so are 1/A, 
A"^ and l/X*). The eigenmode from the top edge 
of the unstable acoustic band was used as an initial 
guess, and the continuation was performed for fre- 
quencies in the gap between the acoustic and optical 
bands. While this led to a family of time-periodic so- 
lutions, they were spatially extended and would not 
be classified as breathers (due to the absence of spa- 
tial localization). These had the same profile as the 
eigenmode from which they were spawned, and the 
amplitude increases as the frequency increases; such 
modes can be classified as nonlinear extensions of 
the (linear) normal modes of the system. Since this 
family bifurcates from the unstable equilibrium, we 
will call it the unstable branch (UB). 

Additionally, when a properly scaled eigenmode 
from the top edge of the stable equilibrium's acous- 
tic band is added to the stable equilibrium, another 
family of spatially extended solutions can be iden- 
tified. This family originates at a frequency equal 
to the top of the stable acoustic band (0.5473), and 
continues down in frequency (up in energy) until it 
collides with UB at a subcritical (in the energy vari- 
able) pitchfork bifurcation. Since this family bifur- 



cates from the stable equilibrium, we will call it the 
stable branch (SB). Sample profiles of the stable and 
unstable branches (SB and UB, respectively), as well 
as the breather branch to which we turn to below, 
are displayed in Fig. S) 

In order to identify localized solutions, we needed 
a suitable initial guess. We excited the shell of the 
first mass within a single lattice site in a chain with a 
total of 200 sites, and then integrated it in the equa- 
tions of motion for a sufficient length of time that 
the extended excitations (phonons) were allowed to 
radiate away, resulting in a localized oscillation cen- 
tered at the excited node. This provided us with a 
sufficiently good initial guess for the profile (trun- 
cated to only = 20 sites) and a reasonable ap- 
proximation to the frequency of an exact breather 
solution. This successfully converged to a breather 
with the Newton-Raphson method, which allowed a 
continuation through the gap. There are many inter- 
esting features to this family of breathers, which we 
call the breather branch (BB) (see again members 
of this family in Fig. U) . Figure [3] shows a view of 
the energy of each solution within this BB as a func- 
tion of frequency and illustrates their corresponding 
stability. It is noticeable that, with the exception 
of a strong instability in the vicinity of the optical 
band and a weaker one in the vicinity of the acous- 
tic band, this branch is generally found to be robust, 
and for a wide parametric regime of frequencies it is 
linearly stable. 

The figures [5] and [5] illustrate the bifurcation 
events associated (directly or indirectly) with the 
breather branch. BB emerges out of a supercriti- 
cal (in the energy variable) pitchfork bifurcation of 
the SB, as is indicated in the left panel of Fig. [S] 
[and thus inherits the (in)stability of this branch]. 
On the other end, for high frequencies, a turning 
point emerges (where the change of monotonicity is 
associated with a stability change) near frequency 
uj = 1.2. Further increase of frequency then leads 
to a rapid decrease in energy, and finally the branch 
merges with the UB in a pitchfork bifurcation that 
is subcritical in the frequency variable but supercrit- 
ical in the energy variable. Notice that the Fig. [5] in 
its right panel also illustrates the subcritical (in the 
energy variable) collision event between the UB and 
SB branches. 

Notice also that in order to visualize each of the 
above mentioned bifurcations and to reveal their 
pitchfork (symmetry-breaking) character, we needed 
to select a proper measure of "asymmetry" A. In 
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FIG. 3. Top: the energy of the members of the breather solution branch is shown as a function of its frequency. 
Notice the logarithmic nature of the vertical scale. Bottom: the Floquet multipliers of each breather solution are 
shown as a function of frequency. Purely real Floquet multipliers deviating from the unit circle are shown in black, 
while complex multipliers with nonzero imaginary part deviating from the unit circle are shown in gray. When neither 
of these deviations from unit modulus arises, the breather branch is linearly stable. 



each case, A was chosen based on how the profiles 
changed before and after the bifurcation. For in- 
stance, let us consider the bifurcation between SB 
and BB. The defining difference between SB and BB 
close to the bifurcation is that solutions in SB have 
a core-1 that is constant across lattice sites, while 
solutions in BB have a core-1 that displays the (as- 
sociated with the breather nature) spatial localiza- 
tion. Here, we chose to have A measure the sample 
skewness of the first core in our system of two cores- 
shells, as follows: 

^= (10) 

where ui is the mean value of core-1. The top- 
left panel of Fig. [5] shows a plot of A vs. energy, 
which explicitly illustrates the supercritical pitch- 
fork structure of the bifurcation. In addition we 
see an exchange in stability of each branch that is 
characteristic of the pitchfork bifurcation. Although 
both branches are unstable before and after the bi- 
furcation (due to a Floquet multiplier with modulus 
greater than one), we see that SB loses a real pair of 
Floquet multipliers, while the BB gains one. 

Next, consider the bifurcation between SB and 



UB. The defining difference between SB and UB 
close to the bifurcation is that solutions in UB have 
homogeneous cores, while solutions in SB have both 
cores (one and two) which are displaced from each 
other. In this case, we choose A to be the mean of all 
four sublattices. The top-right panel of Fig. Oshows 
a plot of A vs. energy, which illustrates the subcriti- 
cal pitchfork bifurcation. Once again, although both 
branches are unstable near the bifurcation, we see 
an exchange of a pair of Floquet multipliers that is 
characteristic of a pitchfork bifurcation; UB loses an 
unstable Floquet multiplier inheriting the stability 
(in this particular subspace) of the SB. 

Lastly, consider the bifurcation between BB and 
UB. Here, the defining difference between the two 
branches close to the bifurcation is that solutions in 
UB have a homogeneous shell-2 in its spatial distri- 
bution, while the spatial localization associated with 
the BB branch leads to an inhomogeneous distribu- 
tion of this field in the latter case. Hence, we have 
chosen A to be the value of V2 at the site with max- 
imum amplitude. The top panel of Fig. [5] shows a 
plot of ^ vs. energy, which illustrates the supercrit- 
ical (in the energy variable) pitchfork bifurcation. 
Lastly, we see in the lower panels of Fig. [B] that BB 
inherits an unstable Floquet multiplier from UB. 



5 




O O O O 0,0 O Q O O 
O O ' O O O O ' o o ^ 



2 

-2 



(b) 



20 

•<l>-#-«>-#-«^ 10- 



OOOOO.OOQOO 

2.„ ' ' « * ' 







3»p»o»o»o»o«o»°» . ♦ . ♦ . ♦°«o«o«o«o«o«o«o^o^ 



10 



15 



20 



"XT 



20 







^ ^ 

^0000 <-,Oa 0000*^0 
- • ♦ • ♦ • ♦ • ♦•♦•♦ •♦• ♦• ♦ •¥• ♦• ♦•♦•♦•♦•♦• ♦ • ♦ • 

o o o o o o o o o oO V 




-20 



o, , ,0 







-©i)QI©B2(i©$©j!®9©» j» • ♦^♦©9s»J©S)S)S©I©i©i©i©i)- 

, o , o , ^ 





4Rf) 
2 




0,0 " ' ' ' 
10 15 
n 



20 



FIG. 4. Sample solutions from each continuation, (a) An extended solution from the stable branch (SB), (b) An 
extended solution from the unstable branch (UB). (c) A nanopteron (phonobreather) from the breather branch (BB). 
This particular solution has a frequency of 0.5989 and a second harmonic which is very close to the second frequency 
in the optical band, (d) A breather solution from BB. This solution has a frequency (1.1936), which is very close 
to the bottom edge of the optical band. Because of this, there are strong resonances with the corresponding optical 
mode, (e) Another breather solution from BB. This solution is stable with a frequency of 1.1280. (f) A third breather 
solution from BB. This solution has a frequency of 0.5884. In these panels, core/shell-1 are shown in blue circles, 
while core/shell-2 are shown as red diamonds. Cores are the filled markers, while shells are open markers. 



As one continues along BB from the lower bifur- 
cation, the energy increases slightly, and then de- 
creases to a local minimum, which occurs very near 
the frequency corresponding to the upper edge of 
the stable acoustic band. Further along BB, we can 
observe in Fig. |3] a series of spikes in the energy of 
the BB branch. These spikes correspond to frequen- 
cies for which the second harmonic of the breather 
is a frequency in the discrete linear band (and are 
thus associated with the corresponding resonances). 
The type of solution that exists around these spikes 
has been termed a phonobreather [16] (also occasion- 
ally referred to as nanopteron due to the resonantly 
excited tail with non- vanishing amplitude). These 
solutions comprise a localized breather solution os- 
cillating at the fundamental frequency, superposed 
with a phonon solution oscillating at the second har- 
monic frequency. 

A sample profile of one such solution is shown 
in Fig. m Since the linear spectrum is so sparse 
for a system of only 20 lattice sites, we can con- 
tinue past these points of resonance by using a step 



size large enough to skip over the exact resonance 
points, and to find breathers for frequencies between 
them. These breathers, termed phantom breathers, 
exist in frequency intervals that should be forbid- 
den for localized solutions in an infinite chain (where 
the linear spectrum becomes dense), as discussed 
more thoroughly in [l6| . Once past the second 
harmonic resonances, there is a long range of fre- 
quencies for which there exist stable, strongly lo- 
calized (breather) solutions. As the frequencies ap- 
proach the bottom edge of the optical band, the en- 
ergy of the solution increases rapidly until it eventu- 
ally reaches the turning point discussed above (and 
changes its stability, as expected by the energy- 
frequency monotonicity criterion; see also 8]) and 
the BB eventually merges with UB as illustrated 
above. 

Having explored the existence and stability prop- 
erties of the discrete breather solutions emerging in 
the nonlinear polarizability model of cores and shells 
in a ferroelectric material, our intention is to comple- 
ment these findings by an examination of the nonlin- 
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FIG. 5. Left: the supercritical (in the energy variable) 
pitchfork bifurcation between SB and BB. The top panel 
shows the relevant 'symmetry' parameter as a function 
of energy. The middle panel shows the Floquet multipli- 
ers of SB as a function of energy, and the bottom panel 
shows the Floquet multipliers of BB as a function of en- 
ergy. Right: the subcritical (again in the energy variable) 
pitchfork bifurcation between UB and SB. The top panel 
shows the 'symmetry' parameter as a function of energy. 
The middle panel shows the Floquet multipliers of SB 
as a function of energy, and the bottom panel shows the 
Floquet multipliers of UB as a function of energy. 
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FIG. 6. The subcritical pitchfork bifurcation between 
BB and UB (for high energies and frequencies near the 
termination point of the breather branch). Top: the 
'symmetry' parameter as a function of energy. Middle: 
the Floquet multipliers of UB as a function of energy. 
Bottom: the Floquet multipliers of BB as a function of 
energy. 



that in a large fraction of the linear gap (and away 
from the edges of their existence region) the BB so- 
lutions are found to be very robust in the long-time 
numerical simulations that we have performed. 



CONCLUSION AND FURTHER WORK 



ear dynamics of such solutions, in the regimes where 
they are found to be unstable. This is particularly 
relevant to explore in the two regimes where strongly 
unstable dynamics emerge, namely for low frequen- 
cies and energies, in the vicinity of the acoustic band, 
as well as for high frequencies and energies near the 
optical band. This is summarized in Fig. [71 In- 
terestingly, it is found that the breather dynamics is 
quite similar in its nature in both of these highly un- 
stable regimes. In particular, we observe that both 
breathers undergo about 25 full oscillations before 
eventually dispersing into linear waves, i.e., the BB 
solutions are completely destroyed by the dynamical 
instability. Nevertheless, we emphasize once again 



In the present work, we have been able to pro- 
vide a direct numerical verification of the existence 
of discrete breathers in the nonlinear model for po- 
larizability that has been proposed for displacive- 
type ferroelectric materials We have identi- 

fied the linear spectrum of the underlying model 
and have shown how breather type excitations can 
emerge within the gaps of the linear spectrum. Ad- 
ditionally, we have provided a detailed continuation 
of such solutions within the gap that revealed a num- 
ber of interesting characteristics. In particular, we 
observed the existence of pitchfork-type symmetry 
breaking bifurcations near the edges of the acous- 
tic and optical bands of the spectrum, that led to 
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FIG. 7. Left: the dynamics of an unstable breather with a frequency of 0.5292. Right: the dynamics of an unstable 
breather with a frequency of 1.1903. In these panels, core/shell-1 are shown in blue circles, while core/shell-2 are 
shown as red diamonds. Cores are the filled markers, while shells are open markers. In each case, the top two panels 
show the profile of the breather on the left and the Floquet multipliers on the right, while the bottom panel shows a 
spacetime plot of the displacements of this breather with a small perturbation. 



the emergence/termination of the breathers, as well 
as the formation of phantom breather solutions (l6| 
in the vicinity of the second harmonic resonances 
of the breather frequency with the linear optical 
band modes. It is also noteworthy that although 
the breathers appear to be destroyed by strong in- 
stabilities associated with real Floquet multipliers in 
the vicinity of the band edges of the gap, they in fact 
appear to be rather stable within a wide parametric 
(frequency) interval in the gap. 

There are many interesting directions along which 
the present study can be extended. On the one hand. 



the two regimes l3, 14 1. In particular, the presence 
of discrete breathers should be reflected in the tem- 
perature and frequency dep endence of vibrational 
and dielectric properties 
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